%% fig3_stability.m
% 本脚本在 I_{ext} 的扫描区间内，
% 利用多个初值求解平衡方程（平衡点写作 [xi, xi, eta]），
% 并利用雅可比矩阵判断各平衡点的稳定性：
%    SNF（稳定节点焦点） —— 用蓝色显示；
%    USF（不稳定鞍焦点） —— 用绿色显示；
%    USP（不稳定鞍点）   —— 用粉显示。
%
% 最后绘制 I_{ext} 与 xi 的散点图，不同稳定性用不同颜色标注，
% 以复刻论文 Fig.3 的平衡轨迹及稳定性分布。

clc; clear; close all;

%% ---- 模型参数设定（与论文一致）----
a      = 0.7;
b      = 0.5;
alpha  = 2.1;
beta   = 2;
gamma  = 2;

%% ---- I_ext 扫描设置 ----
Iext_min = -0.2;
Iext_max =  0.2;
N_Iext   = 201;  % 为获得更细致分叉结构，步长较小
Iext_vals = linspace(Iext_min, Iext_max, N_Iext);

%% 多组初始猜测，每行一组 [xi; eta]
% 这里给出较为常见的初值，尽可能覆盖平衡解可能的初始域
init_guesses = [
    % ======= 偏蓝色分支附近猜测 (xi > 0) =======
     0.1,   0;       % 主体猜测
     0.12,  0.05;    % 小扰动

    % ======= 偏绿色分支附近猜测 (xi < 0) =======
    -0.5,   1.0;
    -0.55,  1.1; 
    -0.6,   1.2;

    % ======= 再多一些 xi 取值略小或更负，eta 在 [0.5, 1.5] 附近 =======
    -0.7,   -1.0;
    -0.2,   0.8;
    -0.8,   1.2;

    % ======= 其他补充，用于捕捉更低分支 (eta < 0) =======
     0.3,  -1.0;
];
nGuess = size(init_guesses,1);

%% ---- 预分配存储 ----
% data_list 中每一行记录: [I_ext, xi, stype]
data_list = [];

%% ---- fsolve 选项 ----
opts = optimset('Display','off','TolFun',1e-12,'TolX',1e-12);

%% ---- 对每个 I_ext 值进行多初值求解，并判断稳定性 ----
for i = 1:N_Iext
    Iext = Iext_vals(i);
    sols = [];      % 用于存储该 I_ext 下的解，每行：[xi, eta]
    types = {};     % 对应的稳定性类型
    for g = 1:nGuess
        guess = init_guesses(g,:)';
        [sol, ~, exitflag] = fsolve(@(X) equilibriumEquations(X, Iext, a, b, alpha, beta, gamma), guess, opts);
        if exitflag > 0
            % 检查重复：判断是否已有解与当前解距离小于 tol
            tol = 1e-3;
            duplicate = false;
            for k = 1:size(sols,1)
                if norm(sols(k,:) - sol') < tol
                    duplicate = true;
                    break;
                end
            end
            if ~duplicate
                sols = [sols; sol']; %#ok<AGROW>
                % 稳定性判断（调用局部函数）
                [~, stype] = checkStabilityEquilibrium(sol, a, b, alpha, beta, gamma);
                types{end+1} = stype;
            end
        end
    end
    % 如果找到解，则保存每个解对应的数据
    if ~isempty(sols)
        nSol = size(sols,1);
        for j = 1:nSol
            data_list = [data_list; Iext, sols(j,1), type2num(types{j})];
        end
    end
end

%% 分离数据：以不同颜色绘图
% data_list 列分别为 I_ext, xi, 类型编号
I_data = data_list(:,1);
xi_data = data_list(:,2);
stype_num = data_list(:,3);

% 定义颜色映射：SNF - 蓝色, USF - 绿色, USP - 洋红色, 其他 - 黑色
colors = {'b','g','m','k'};
% 对应稳定性编号：用辅助函数 type2num，返回 1: SNF, 2: USF, 3: USP; 若其他返回 4
figure;
hold on; box on;
for j = 1:length(I_data)
    % 绘制每个点
    color = colors{stype_num(j)};
    plot(I_data(j), xi_data(j), '.', 'Color', color, 'MarkerSize', 16);
end
xlabel('I_{ext}','FontSize',12);
ylabel('\xi','FontSize',12);
title('Fig.3 平衡轨迹及稳定性','FontSize',14);
grid on;
xlim([Iext_min Iext_max]);
legend({'SNF (蓝)','USF (绿)','USP (粉)'}, 'Location','best');

%% ---- 局部函数：平衡方程 ----
function F = equilibriumEquations(X, Iext, a, b, alpha, beta, gamma)
    % 对于平衡点 s = [xi, xi, eta]，令
    % Vx = (alpha+beta*xi)*xi,
    % 平衡条件为：
    %   F1 = gamma*(alpha+beta*xi)*xi + tanh(eta) - eta = 0,
    %   F2 = Iext - (eta^2 - 0.5)*Vx + a*abs(Vx) + b*Vx = 0.
    xi  = X(1);
    eta = X(2);
    Vx = (alpha + beta*xi)*xi;
    F1 = gamma*(alpha+beta*xi)*xi + tanh(eta) - eta;
    F2 = Iext - (eta^2 - 0.5)*Vx + a*abs(Vx) + b*Vx;
    F = [F1; F2];
end

%% ---- 局部函数：稳定性判断 ----
function [lambda_vals, stype] = checkStabilityEquilibrium(X, a, b, alpha, beta, gamma)
    % X 为解向量 [xi; eta]，平衡点写为 [xi, xi, eta]
    xi  = X(1);
    eta = X(2);
    Vx = (alpha+beta*xi)*xi;
    GM = eta^2 - 0.5;
    % 对于 f1: 其对 q 的偏导（q=xi），注意 q = sigma = xi
    dVx_dq = alpha + beta*xi;
    dVx_dsigma = beta*xi;
    % 计算 f1 关于 q, sigma, phi 的偏导数
    J11 = -GM * dVx_dq + (a*sign(Vx)+b)*dVx_dq;
    J12 = -GM * dVx_dsigma + (a*sign(Vx)+b)*dVx_dsigma;
    J13 = -2*eta*Vx;
    % f2 = q - sigma, 对 q, sigma, phi 的偏导：J21=1, J22=-1, J23=0
    J21 = 1; J22 = -1; J23 = 0;
    % f3 = gamma*Vx + tanh(eta) - eta，对 q, sigma, phi 的偏导：
    J31 = gamma*dVx_dq;
    J32 = gamma*dVx_dsigma;
    J33 = sech(eta)^2 - 1;
    J = [J11, J12, J13;
         J21, J22, J23;
         J31, J32, J33];
    lambda_vals = eig(J);
    % 根据特征值判断稳定性：若所有实部均 < 0 则 SNF；若恰有 1 个正则 USF；否则 USP
    rp = real(lambda_vals);
    if all(rp < 0)
        stype = 'SNF';
    elseif (sum(rp > 0) == 1 && sum(rp < 0) == 2)
        stype = 'USF';
    else
        stype = 'USP';
    end
end

%% ---- 局部函数：稳定性类型转换 ----
function num = type2num(stype)
    % 将稳定性类型 stype 转换为数字编号，用于颜色映射：
    % 'SNF' -> 1, 'USF' -> 2, 'USP' -> 3, 其他 -> 4
    if strcmp(stype, 'SNF')
        num = 1;
    elseif strcmp(stype, 'USF')
        num = 2;
    elseif strcmp(stype, 'USP')
        num = 3;
    else
        num = 4;
    end
end
